CapStone Project_2025
Slides: slides.html ( Go to slides.qmd to edit)
Diabetes mellitus (DM) represents a major public health concern and is closely associated with factors such as obesity, age, race, and gender. Identifying these associated risk factors is essential for developing targeted and effective intervention strategies D’Angelo and Ran (2025). Logistic Regression that estimates the association between risk factors and binary outcomes (presence or absence of diabetes) is a standard analytical approach but, is insufficient in analyzing the complexity of healthcare data (DNA sequences, imaging, patient-reported outcomes, electronic health records (EHRs), longitudinal health measurements, diagnoses, and treatments. Zeger et al. (2020). Classical maximum likelihood estimation (MLE) yields unstable results in small samples with missing data, or quasi- and complete separation.
The Bayesian hierarchical model with Markov Chain Monte Carlo (MCMC) on multivariate longitudinal healthcare data with two levels: (1) repeated measures over time within individuals and (2) individuals nested within a population, by integrating prior knowledge and by added exogenous covariates (e.g., age, clinical history), and endogenous covariates (e.g., current treatment), yield posterior distributions and marginal distributions from MCMC estimation of parameters provide risk prediction (pneumonia, prostate cancer, and mental disorders). The model’s limitation is its parametric nature.
In Bayesian inference Chatzimichail and Hatjimihail (2023), comparisons between parametric models (with fixed sets of parameters) and non-parametric models (which make no prior assumptions about the distribution) applied to the National Health and Nutrition Examination Survey demonstrated the use of posterior probabilities to classify diseases. In contrast, conventional approaches that rely on fixed clinical thresholds and dichotomous classifications often fail to capture the complex relationships between diagnostic tests and disease prevalence. These traditional methods tend to oversimplify the underlying heterogeneity across populations and perform poorly when dealing with skewed, bimodal, or multimodal data distributions.
Bayesian nonparametric modeling is flexible, adaptable, and robust alternative to traditional parametric approaches capturing complex and multimodal data patterns by integration of priors with multiple diagnostic tests enhances diagnostic accuracy and precision. However, the approach can be limited by restricted access to scholarly resources and excessive reliance on prior information. Combined with other statistical and computational techniques, Bayesian nonparametric models can further strengthen diagnostic performance addressing systemic bias, incomplete data, and non-representative or non-normal datasets Chatzimichail and Hatjimihail (2023)
The Bayesian framework described by Schoot et al. (2021) highlights the critical role of priors, data modeling, inference, model checking, posterior sampling, variational inference, and variable selection in applications across social sciences, ecology, genetics, and medicine. Variable selection is important because multicollinearity, limited sampling, and overfitting can lead to weak predictive performance and reduced interpretability. Based on the degree of uncertainty in hyperparameters, priors are classified as informative, weakly informative, or diffuse (larger variances indicating greater uncertainty).
Prior elicitation may be based on expert opinion, generic knowledge, or data-driven using maximum likelihood estimates or sample statistics. Prior sensitivity analysis assesses the alignment between priors and the likelihood. Prior improve estimation via data-informed shrinkage, regularization, and influence on algorithms by focusing on high-density regions. Specifying prior information strengthens inference by integrating observed data with uncertain parameters in small or less informative samples, generating a range of possible values through the likelihood function. The Markov Chain Monte Carlo (MCMC) simulation algorithm, ( brms and blavaan in R), allows empirical estimation of posterior distributions. Spatial and temporal Bayesian models analyze large-scale cancer genomics, identify novel molecular interactions, detection of mutational signatures, genomic-based patient stratification for clinical trials, and insights into cancer evolution, but ahve limited application because of temporal autocorrelation and the inherent subjectivity in prior elicitation.
Bayesian normal linear regression, has been applied in metrology for instrument calibration using parametric (conjugate) models with a Normal–Inverse-Gamma (NIG) prior Klauenberg et al. (2015) as in Gaussian regression, errors are independent and identically distributed with unknown variance, the relationship between X and Y is statistical and due to the associated noise and model uncertainty, can not be treated as a direct measurement function. The Bayesian method accounts for uncertainty, treats both observables (data) and unobservables (parameters and auxiliary variables) as random and assigns probability distributions to unobservables. Incorporation of prior, enhances interpretability and robustness. The NIG distribution (conjugate prior) enables the use of vague or non-informative priors, and hierarchical priors further add flexibility by modeling uncertainty across multiple levels, and capture complex data nuances.
Bayesian Hierarchical / meta-analytic linear regression model augments data by incorporating both exchangeable and unexchangeable information on parameters, addressing issues associated with multiple testing providing low statistical power when separate significance tests are conducted across studies with different predictors. Adding priors from meta-analysis in Bayesian Hierarchical regression address the challenge of small sample size with limited availability of previous articles, and resolve limitation of univariate analyses, and the relationship issues among multiple regression parameters within a study. (1) Exchangeable Priors are when the current data and previous studies share the same set of predictors, and (2) Unexchangeable when the predictors differ. The hierarchical unexchangeable model is applicable in studying differences in studies, enabling explicit testing of the exchangeability assumption. Approach is limited due to the correlation between identical set of predictors. Leeuw and Klugkist (2012)
A sequential clinical reasoning model** Liu et al. (2013) demonstrated screening adults (20–79 years, Taiwan) while addressing the limited availability of molecular information and provided an alternative method leveraging routinely collected biological markers classifying diseases. Sequential adding of predictors in three models: (1) demographic features (basic model), (2) six metabolic components (metabolic score model), and (3) conventional risk factors (enhanced model), incorporating priors, and emulating a clinician’s evaluation process, the model assumes normally distributed regression coefficients, accounts for uncertainty in clinical weights, and averages credible intervals for predicted risk estimates. by incorporating patient background, individual characteristics, in the enhanced model, significantly contributed to the baseline risk estimation by capturing ecological heterogeneity. The model limitation is interactions between predictors and the external cross-validation issue.
Bayesian multiple imputation with logistic regression addresses missing data in clinical research Austin et al. (2021) by classifying and analyzing missing values based on (i) patients’ refusal (ii) loss to follow-up, (iii) investigator or mechanical errors, or (iv) physicians’ choice on ordering investigations. Missing at random (MAR), missing not at random (MNAR), or missing completely at random (MCAR) are crucial in data analysis. Multiple imputation (MI) (R, SAS, or Stata) provides plausible values while simultaneously conducting identical statistical analyses across them, creating multiple completed datasets, and the pooled results classify patients’ health status and mortality rates.
Aims
The study aims to apply Bayesian logistic regression to predict diabetes status and to evaluate the associations between body mass index (BMI), age (≥20 years), gender, and race as predictors, using a retrospective dataset from the 2013–2014 NHANES survey. NHANES employs a complex sampling design, including stratification, clustering, and oversampling of specific population subgroups, rather than uniform random sampling. A Bayesian analytical approach is used to address challenges posed by dataset anomalies such as missing data, complete case analysis, and separation that limit the efficiency and reliability of traditional logistic regression in predicting health outcomes.
Statistical Tool: we use R, R packages and libraries to import data, perform data wrangling and analysis.
Data source: NHANES 2-year data is a cross-sectional weighted data (2013-2014 year) Center for Health Statistics (1999). Three datasets (demographics, exam, questionnaire) imported (Haven package to coverted .XPT files in R to dataframe (df)) and after selecting variables of interest a merged dataset was created.
Modeling - We conduct bayesian logisitc regression to estimate the association between BMI, age, sex, and race/ethnicity and predict doctor-diagnosed diabetes (DIQ010).
# loading packages
options(repos = c(CRAN = "https://cloud.r-project.org"))
install.packages("nhanesA")
library ("nhanesA")
library(tidyverse)
library(knitr)
library(ggthemes)
library(ggrepel)
library(dslabs)
library(Hmisc)
library(dplyr)
library(tidyr)
library(forcats)
library(ggplot2)
library(classpackage)
library(janitor)
install.packages("gt")
library(gt)
library(survey)
library(DataExplorer)
library(logistf)Method
Data pre-processing and cleaning Subsets are created from the original weighted datasets (demographics, exam, questionnaire) with selected variables are merged using ID to create a single dataframe.
DIQ050 a secondary variable describing treatment status (insulin use) to exclude those casesMerged dataset created, cleaned and visualized for any anomalies and ata Exploration.
# weighted means of each variable
str(merged_data)'data.frame': 10175 obs. of 10 variables:
$ SEQN : num 73557 73558 73559 73560 73561 ...
$ RIDAGEYR: num 69 54 72 9 73 56 0 61 42 56 ...
$ RIAGENDR: Factor w/ 2 levels "Male","Female": 1 1 1 1 2 1 1 2 1 2 ...
$ RIDRETH1: Factor w/ 5 levels "Mexican American",..: 4 3 3 3 3 1 3 3 2 3 ...
$ SDMVPSU : num 1 1 1 2 2 1 1 1 2 1 ...
$ SDMVSTRA: num 112 108 109 109 116 111 105 114 106 112 ...
$ WTMEC2YR: num 13481 24472 57193 55767 65542 ...
$ BMXBMI : num 26.7 28.6 28.9 17.1 19.7 41.7 NA 35.7 NA 26.5 ...
$ DIQ010 : Factor w/ 5 levels "Yes","No","Borderline",..: 1 1 1 2 2 2 NA 2 2 2 ...
$ DIQ050 : Factor w/ 4 levels "Yes","No","Refused",..: 1 1 1 2 2 2 NA 2 2 2 ...plot_str(merged_data)
introduce(merged_data) rows columns discrete_columns continuous_columns all_missing_columns
1 10175 10 4 6 0
total_missing_values complete_rows total_observations memory_usage
1 1933 9053 101750 656968p1 <- plot_intro(merged_data, title="Figure 1 (Merged dataset). Structure of variables and missing observations.")
p2 <- plot_missing(merged_data, title="Figure 2(Merged dataset). Breakdown of missing observations.")
# Save it as a PNG file
ggsave("Figure1_MergedDataset.png", plot = p1, width = 8, height = 6, dpi = 300)
ggsave("Figure2_MergedDataset.png", plot = p2, width = 8, height = 6, dpi = 300)# print(glimpse(merged_data))
print(table(merged_data$BMDBMIC, useNA = "ifany"))< table of extent 0 >print(table(merged_data$DIQ010, useNA = "ifany"))
Yes No Borderline Refused Don't know <NA>
737 8841 185 1 5 406 # ---- Coercion helpers (handle labelled/character) ----
to_num <- function(x) {
if (is.numeric(x)) return(x)
xc <- as.character(x)
n <- suppressWarnings(readr::parse_number(xc))
if (mean(is.na(n)) > 0.80) {
xlow <- tolower(trimws(xc))
n <- dplyr::case_when(
xlow %in% c("1","yes","yes, told") ~ 1,
xlow %in% c("2","no","no, not told") ~ 2,
xlow %in% c("3","borderline") ~ 3,
xlow %in% c("7","refused") ~ 7,
xlow %in% c("9","don't know","dont know","unknown") ~ 9,
TRUE ~ NA_real_
)
}
as.numeric(n)
}
merged_data <- merged_data %>%
mutate(
DIQ010 = to_num(DIQ010),
DIQ050 = to_num(if (!"DIQ050" %in% names(.)) NA_real_ else DIQ050),
BMXBMI = suppressWarnings(as.numeric(BMXBMI)),
RIDAGEYR = suppressWarnings(as.numeric(RIDAGEYR)),
RIAGENDR = suppressWarnings(as.numeric(RIAGENDR)),
RIDRETH1 = suppressWarnings(as.numeric(RIDRETH1)),
SDMVPSU = suppressWarnings(as.numeric(SDMVPSU)),
SDMVSTRA = suppressWarnings(as.numeric(SDMVSTRA)),
WTMEC2YR = suppressWarnings(as.numeric(WTMEC2YR))
)
# ---- Diagnostics BEFORE save ----
cat("DIQ010 counts BEFORE save:\n")DIQ010 counts BEFORE save:print(table(merged_data$DIQ010, useNA = "ifany"))
1 2 3 7 9 <NA>
737 8841 185 1 5 406 cat("Count with DIQ010 in {1,2}:", sum(merged_data$DIQ010 %in% c(1,2), na.rm = TRUE), "\n")Count with DIQ010 in {1,2}: 9578 # ---- Save to file for reuse ----
dir.create("data", showWarnings = FALSE)
# ---- Save ----
dir.create("data", showWarnings = FALSE, recursive = TRUE)
saveRDS(merged_data, "data/merged_2013_2014.rds")
message("Saved: data/merged_2013_2014.rds")EDA
##
# ---------------- Basic Exploration (adults) ----------------
# Keep adults only and build analysis variables
adult <- merged_data %>%
dplyr::filter(RIDAGEYR >= 20) %>%
dplyr::transmute(
# --- keep survey design variables so svydesign() can see them ---
SDMVPSU, SDMVSTRA, WTMEC2YR,
# --- outcome: DIQ010 (1 yes, 2 no; 3/7/9 -> NA) ---
diabetes_dx = dplyr::case_when(
DIQ010 == 1 ~ 1,
DIQ010 == 2 ~ 0,
DIQ010 %in% c(3, 7, 9) ~ NA_real_,
TRUE ~ NA_real_
),
# --- predictors (raw) ---
bmi = BMXBMI,
age = RIDAGEYR,
# sex (1=Male, 2=Female)
sex = forcats::fct_recode(factor(RIAGENDR), Male = "1", Female = "2"),
# race (5-level)
race = forcats::fct_recode(
factor(RIDRETH1),
"Mexican American" = "1",
"Other Hispanic" = "2",
"NH White" = "3",
"NH Black" = "4",
"Other/Multi" = "5"
),
# keep DIQ050 so we can safely reference it (may be absent/NA in some rows)
DIQ050 = DIQ050
) %>%
# standardize continuous predictors
dplyr::mutate(
age_c = as.numeric(scale(age)),
bmi_c = as.numeric(scale(bmi)),
bmi_cat = cut(
bmi,
breaks = c(-Inf, 18.5, 25, 30, 35, 40, Inf),
labels = c("<18.5","18.5–<25","25–<30","30–<35","35–<40","≥40"),
right = FALSE
)
) %>%
# adjust outcome: if female & DIQ050==1 ("only when pregnant"), set to 0 (not diabetes)
dplyr::mutate(
diabetes_dx = ifelse(sex == "Female" & !is.na(DIQ050) & DIQ050 == 1, 0, diabetes_dx)
)
# Make NH White the reference level for race (clearer interpretation)
adult <- adult %>%
dplyr::mutate(
race = forcats::fct_relevel(race, "NH White")
)
# --- sanity checks ---
cat("Adults n =", nrow(adult), "\n")Adults n = 5769 # data exploration
print(table(adult$diabetes_dx, useNA = "ifany"))
0 1 <NA>
4974 618 177 print(table(adult$sex, useNA = "ifany"))
Male Female
2758 3011 print(table(adult$race, useNA = "ifany"))
NH White Mexican American Other Hispanic NH Black
2472 767 508 1177
Other/Multi
845 if (sum(!is.na(adult$diabetes_dx)) == 0) {
stop("Too few non-missing outcomes for modeling (n = 0). Check DIQ010 upstream.")
}
# (optional plots omitted for brevity)
# save for downstream
if (!dir.exists("data")) dir.create("data", recursive = TRUE)
saveRDS(adult, "data/adult_cleaned_2013_2014.rds")# survey design
# ---------------- Survey Design ----------------
# Use exam weights because BMI (BMXBMI) is an MEC variable
nhanes_design_adult <- survey::svydesign(
id = ~SDMVPSU,
strata = ~SDMVSTRA,
weights = ~WTMEC2YR,
nest = TRUE,
data = adult
)
# quick weighted checks
survey::svymean(~age, nhanes_design_adult, na.rm = TRUE) mean SE
age 47.496 0.3805survey::svymean(~diabetes_dx, nhanes_design_adult, na.rm = TRUE) mean SE
diabetes_dx 0.089016 0.0048** Statistical Methods**
Multiple Logistic regression on survey weighted dataset
-We conducted frequentist method Multiple Logistic regression on a survey-weighted dataset, for complete case analysis and data exploration
Multivariate Imputation by Chained Equations (MICE) - We then conducted MICE to manage missiging data - Considering the small sample size, Multivariate Imputation by Chained Equations (MICE) was conducted as an alternative to the Bayesian Approach Buuren and Groothuis-Oudshoorn (2011)
Multiple imputation (MI) is an alternative analytic approach for small dataset with missingness.
Flatness of the density, heavy tails, non-zero peakedness, skewness and multimodality do not hamper the good performance of multiple imputation for the mean structure in samples n > 400 even for high percentages (75%) of missing data in one variable @Van Buuren and Van Buuren (2012).
Multiple Imputation (MI) can be performed using mice package in R and adds sampling variability to the imputations.
Iterative MICE imputes missing values of one variable at a time, using regression models based on the other variables in the dataset.
In the chain process, each imputed variable become a predictor for the subsequent imputation, and the entire process is repeated multiple times to create several complete datasets, each reflecting different possibilities for the missing data.
Bayesian Logistic Regression
Bayesian statistics is about updating beliefs with evidence:
Posterior ∝ Likelihood × Prior
Prior (p(θ)): Your initial belief about a parameter before seeing the data.
Likelihood (p(y|θ)): How probable the observed data are given the parameters. This is derived from the model (e.g., logistic regression likelihood).
Posterior (p(θ|y)): Your updated belief about the parameter after seeing the data.
We selected Bayesian Logistic Regression since our study response variable is a binary outcome (Diabetes:yes/no)
Bayesian Logistic Regression is based on binomial probability Bayes’ rules, and predicts probability of disease outcome
Bayes analyzes linear relation between the predictor (Age, Race, BMI, Gender) and outcome response variable (Diabetes).
it considers that predictors and response variables are independent.
Regression a of a discrete variable (0 or 1) is a Bernoulli probability model that classifies categorical response variables - predicting Diabetes.
Logit link provides probabilities for the response variable.
We use Weakly informative priors Normal (0, 2.5) for logistic regression coefficients and intercept, Normal(0, 10), allows a wide range of baseline log-odds and helps with convergence and avoids extreme estimates. Good default for most applications in social, health, or epidemiological studies.
In Bayesian statistics, every unknown parameter (like a regression coefficient, mean, or variance) is treated as a random variable with a probability distribution that reflects uncertainty.
we then summarize in Bayesian results (posterior mean, credible intervals, etc.).
Bayesian Logistic Regression model
\[ \text{logit}(P(Y_i=1)) = \beta_0 + \beta_1 \cdot Age_i + \beta_2 \cdot BMI_i + \beta_3 \cdot Race_i + \beta_4 \cdot Gender_i \]
Linear Regression equation:
\[ y_i = \beta_0 + \beta_1 X_1 +\varepsilon_i \]
Diagnostics performed before Bayesian Regression
# Modeling
library(broom)
library(mice)
library(brms)
library(posterior)
library(bayesplot)
library(knitr)
# --- Guardrails for modeling ---
n_outcome <- sum(!is.na(adult$diabetes_dx))
if (n_outcome == 0) stop("Too few non-missing outcomes for modeling. n = 0")
# Ensure factors and >=2 observed levels among complete outcomes
adult <- adult %>%
dplyr::mutate(
sex = if (!is.factor(sex)) factor(sex) else sex,
race = if (!is.factor(race)) factor(race) else race
)
if (nlevels(droplevels(adult$sex[!is.na(adult$diabetes_dx)])) < 2)
stop("sex has <2 observed levels after filtering; check data availability.")
if (nlevels(droplevels(adult$race[!is.na(adult$diabetes_dx)])) < 2)
stop("race has <2 observed levels after filtering; check Data Prep.")
# ------------------------- 1) Survey-weighted complete-case -------------------------
# Build a logical filter on the original adult data (same length as design$data)
keep_cc <- with(
adult,
!is.na(diabetes_dx) & !is.na(age_c) & !is.na(bmi_c) &
!is.na(sex) & !is.na(race)
)
# Subset the survey design using the logical vector (same length as original)
des_cc <- subset(nhanes_design_adult, keep_cc)
# Corresponding complete-case data (optional)
cc <- adult[keep_cc, ] |> droplevels()
cat("\nComplete-case N for survey-weighted model:", nrow(cc), "\n")
Complete-case N for survey-weighted model: 5349 print(table(cc$race))
NH White Mexican American Other Hispanic NH Black
2293 713 470 1101
Other/Multi
772 print(table(cc$diabetes_dx))
0 1
4752 597 print(table(cc$sex))
Male Female
2551 2798 form_cc <- diabetes_dx ~ age_c + bmi_c + sex + race
svy_fit <- survey::svyglm(formula = form_cc, design = des_cc, family = quasibinomial())
summary(svy_fit)
Call:
svyglm(formula = form_cc, design = des_cc, family = quasibinomial())
Survey design:
subset(nhanes_design_adult, keep_cc)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.67143 0.11935 -22.383 1.68e-08 ***
age_c 1.10833 0.05042 21.981 1.94e-08 ***
bmi_c 0.63412 0.05713 11.099 3.88e-06 ***
sexFemale -0.63844 0.10926 -5.843 0.000386 ***
raceMexican American 0.71091 0.13681 5.196 0.000826 ***
raceOther Hispanic 0.46469 0.13474 3.449 0.008712 **
raceNH Black 0.51221 0.15754 3.251 0.011677 *
raceOther/Multi 0.84460 0.17756 4.757 0.001433 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasibinomial family taken to be 0.8455444)
Number of Fisher Scoring iterations: 6plot(residuals(svy_fit, type='deviance'))
# Survey-weighted OR table (no intercept)
svy_or <- broom::tidy(svy_fit, conf.int = TRUE) %>%
dplyr::mutate(OR = exp(estimate), LCL = exp(conf.low), UCL = exp(conf.high)) %>%
dplyr::select(term, OR, LCL, UCL, p.value) %>%
dplyr::filter(term != "(Intercept)")
knitr::kable(svy_or, caption = "Survey-weighted odds ratios (per 1 SD)")| term | OR | LCL | UCL | p.value |
|---|---|---|---|---|
| age_c | 3.0292807 | 2.6967690 | 3.4027912 | 0.0000000 |
| bmi_c | 1.8853571 | 1.6526296 | 2.1508579 | 0.0000039 |
| sexFemale | 0.5281132 | 0.4104905 | 0.6794397 | 0.0003857 |
| raceMexican American | 2.0358434 | 1.4850041 | 2.7910081 | 0.0008262 |
| raceOther Hispanic | 1.5915182 | 1.1664529 | 2.1714810 | 0.0087119 |
| raceNH Black | 1.6689718 | 1.1605895 | 2.4000450 | 0.0116773 |
| raceOther/Multi | 2.3270527 | 1.5451752 | 3.5045697 | 0.0014331 |
The residual plot looks okay and does not show any pattern.
fit_mi <- with(imp, {
age_c <- as.numeric(scale(age))
bmi_c <- as.numeric(scale(bmi))
glm(diabetes_dx ~ age_c + bmi_c + sex + race, family = binomial())
})
pool_mi <- pool(fit_mi)
summary(pool_mi) term estimate std.error statistic df p.value
1 (Intercept) -2.6895645 0.09941301 -27.054453 5566.204 1.486581e-151
2 age_c 1.0660265 0.05594733 19.054108 5520.446 1.911564e-78
3 bmi_c 0.5468538 0.04473386 12.224604 5148.557 6.751227e-34
4 sexFemale -0.6178297 0.09379129 -6.587282 5551.660 4.892566e-11
5 raceMexican American 0.8877355 0.13750463 6.456041 5472.583 1.167455e-10
6 raceOther Hispanic 0.5606621 0.17485537 3.206433 5573.987 1.351505e-03
7 raceNH Black 0.6809629 0.11981185 5.683602 5576.734 1.385727e-08
8 raceOther/Multi 0.7476406 0.15300663 4.886328 4749.963 1.061140e-06## table
mi_or <- summary(pool_mi, conf.int = TRUE, exponentiate = TRUE) %>%
dplyr::rename(
term = term, OR = estimate, LCL = `2.5 %`, UCL = `97.5 %`, p.value = p.value
) %>%
dplyr::filter(term != "(Intercept)")
knitr::kable(mi_or, caption = "MI pooled odds ratios (per 1 SD)")| term | OR | std.error | statistic | df | p.value | LCL | UCL | conf.low | conf.high | |
|---|---|---|---|---|---|---|---|---|---|---|
| 2 | age_c | 2.9038183 | 0.0559473 | 19.054108 | 5520.446 | 0.0000000 | 2.6021752 | 3.2404277 | 2.6021752 | 3.2404277 |
| 3 | bmi_c | 1.7278084 | 0.0447339 | 12.224604 | 5148.557 | 0.0000000 | 1.5827382 | 1.8861754 | 1.5827382 | 1.8861754 |
| 4 | sexFemale | 0.5391132 | 0.0937913 | -6.587282 | 5551.660 | 0.0000000 | 0.4485669 | 0.6479368 | 0.4485669 | 0.6479368 |
| 5 | raceMexican American | 2.4296216 | 0.1375046 | 6.456041 | 5472.583 | 0.0000000 | 1.8555327 | 3.1813298 | 1.8555327 | 3.1813298 |
| 6 | raceOther Hispanic | 1.7518320 | 0.1748554 | 3.206433 | 5573.987 | 0.0013515 | 1.2434346 | 2.4680953 | 1.2434346 | 2.4680953 |
| 7 | raceNH Black | 1.9757793 | 0.1198118 | 5.683602 | 5576.734 | 0.0000000 | 1.5621842 | 2.4988753 | 1.5621842 | 2.4988753 |
| 8 | raceOther/Multi | 2.1120110 | 0.1530066 | 4.886328 | 4749.963 | 0.0000011 | 1.5646727 | 2.8508138 | 1.5646727 | 2.8508138 |
library(gt)
# 3) Bayesian Logistic Regression (formula weights)
adult_imp1 <- complete(imp, 1) %>%
dplyr::mutate(
age_c = as.numeric(scale(age)),
bmi_c = as.numeric(scale(bmi)),
wt_norm = WTMEC2YR / mean(WTMEC2YR, na.rm = TRUE),
# ensure factor refs match survey/MICE:
race = forcats::fct_relevel(race, "NH White"),
sex = forcats::fct_relevel(sex, "Male")
) %>%
dplyr::filter(!is.na(diabetes_dx), !is.na(age_c), !is.na(bmi_c),
!is.na(sex), !is.na(race)) %>%
droplevels()
stopifnot(all(is.finite(adult_imp1$wt_norm)))
priors <- c(
set_prior("normal(0, 2.5)", class = "b"),
set_prior("student_t(3, 0, 10)", class = "Intercept")
)
bayes_fit <- brm(
formula = diabetes_dx | weights(wt_norm) ~ age_c + bmi_c + sex + race,
data = adult_imp1,
family = bernoulli(link = "logit"),
prior = priors,
chains = 4, iter = 2000, seed = 123,
control = list(adapt_delta = 0.95),
refresh = 0 # quiet Stan output
)Running /opt/R/4.4.2/lib/R/bin/R CMD SHLIB foo.c
using C compiler: ‘gcc (GCC) 11.5.0 20240719 (Red Hat 11.5.0-2)’
gcc -I"/opt/R/4.4.2/lib/R/include" -DNDEBUG -I"/opt/R/4.4.2/lib/R/library/Rcpp/include/" -I"/opt/R/4.4.2/lib/R/library/RcppEigen/include/" -I"/opt/R/4.4.2/lib/R/library/RcppEigen/include/unsupported" -I"/opt/R/4.4.2/lib/R/library/BH/include" -I"/opt/R/4.4.2/lib/R/library/StanHeaders/include/src/" -I"/opt/R/4.4.2/lib/R/library/StanHeaders/include/" -I"/opt/R/4.4.2/lib/R/library/RcppParallel/include/" -I"/opt/R/4.4.2/lib/R/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/opt/R/4.4.2/lib/R/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fpic -g -O2 -c foo.c -o foo.o
In file included from /opt/R/4.4.2/lib/R/library/RcppEigen/include/Eigen/Core:19,
from /opt/R/4.4.2/lib/R/library/RcppEigen/include/Eigen/Dense:1,
from /opt/R/4.4.2/lib/R/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
from <command-line>:
/opt/R/4.4.2/lib/R/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
679 | #include <cmath>
| ^~~~~~~
compilation terminated.
make: *** [/opt/R/4.4.2/lib/R/etc/Makeconf:195: foo.o] Error 1summary(bayes_fit) Family: bernoulli
Links: mu = logit
Formula: diabetes_dx | weights(wt_norm) ~ age_c + bmi_c + sex + race
Data: adult_imp1 (Number of observations: 5592)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -2.66 0.09 -2.84 -2.50 1.00 4187 3510
age_c 1.09 0.06 0.97 1.22 1.00 3012 3098
bmi_c 0.63 0.05 0.53 0.72 1.00 3472 3315
sexFemale -0.66 0.10 -0.86 -0.46 1.00 4003 3052
raceMexicanAmerican 0.69 0.18 0.35 1.04 1.00 3526 2843
raceOtherHispanic 0.43 0.24 -0.07 0.89 1.00 4058 3114
raceNHBlack 0.54 0.15 0.24 0.82 1.00 3597 3177
raceOtherDMulti 0.82 0.19 0.45 1.19 1.00 3763 3257
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).prior_summary(bayes_fit) prior class coef group resp dpar nlpar lb ub
normal(0, 2.5) b
normal(0, 2.5) b age_c
normal(0, 2.5) b bmi_c
normal(0, 2.5) b raceMexicanAmerican
normal(0, 2.5) b raceNHBlack
normal(0, 2.5) b raceOtherDMulti
normal(0, 2.5) b raceOtherHispanic
normal(0, 2.5) b sexFemale
student_t(3, 0, 10) Intercept
tag source
user
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
user###
library(ggplot2)
# Example: priors for two coefficients (age and bmi)
prior_draws <- tibble(
term = rep(c("Age (per 1 SD)", "BMI (per 1 SD)"), each = 4000),
value = c(rnorm(4000, 0, 2.5), rnorm(4000, 0, 2.5))
)
ggplot(prior_draws, aes(x = value, fill = term)) +
geom_density(alpha = 0.5) +
theme_minimal() +
labs(title = "Prior Distributions for Coefficients",
x = "Coefficient Value", y = "Density") +
scale_fill_manual(values = c("skyblue", "orange"))
# Diabetes vs BMI
library(ggplot2)
# Create the plot
p3 <- ggplot(adult_imp1, aes(x = factor(diabetes_dx), y = bmi, fill = factor(diabetes_dx))) +
geom_boxplot(alpha = 0.7) +
scale_x_discrete(labels = c("0" = "No Diabetes", "1" = "Diabetes")) +
labs(
x = "Diabetes Diagnosis",
y = "BMI",
title = "BMI Distribution by Diabetes Status"
) +
theme_minimal() +
theme(legend.position = "none")
# Save the plot as an image file (PNG)
ggsave("BMI_Distribution_by_Diabetes_Status.png", plot = p3, width = 7, height = 5, dpi = 300)
p3
# logistic regression curve
p4 <- ggplot(adult_imp1, aes(x = bmi, y = diabetes_dx)) +
geom_point(aes(y = diabetes_dx), alpha = 0.2, position = position_jitter(height = 0.02)) +
geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE, color = "blue") +
labs(
x = "BMI",
y = "Probability of Diabetes",
title = "Predicted Probability of Diabetes vs BMI"
) +
theme_minimal()
ggsave("Predicted Probability of Diabetes vs BMI.png", plot = p4, width = 7, height = 5, dpi = 300)
p4
Once we get posterior draws, we study Summary stats Mean, median, 95% credible intervals summary(bayes_fit) or posterior_summary(bayes_fit) - Visualization Distribution shape mcmc_hist(posterior, pars = c(“b_age”)) or mcmc_areas() - Pairwise plots Correlations between parameters mcmc_pairs(posterior) - Posterior predictive checks Compare model predictions vs observed data pp_check(bayes_fit) - Model comparison Using LOO or WAIC loo(bayes_fit) or waic(bayes_fit)
Assumptions for Bayesian logistic regression - posterior check - plots for linearity - mcmc trace plots for convergence - bayes_R2 for model fit
summary(bayes_fit) Family: bernoulli
Links: mu = logit
Formula: diabetes_dx | weights(wt_norm) ~ age_c + bmi_c + sex + race
Data: adult_imp1 (Number of observations: 5592)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -2.66 0.09 -2.84 -2.50 1.00 4187 3510
age_c 1.09 0.06 0.97 1.22 1.00 3012 3098
bmi_c 0.63 0.05 0.53 0.72 1.00 3472 3315
sexFemale -0.66 0.10 -0.86 -0.46 1.00 4003 3052
raceMexicanAmerican 0.69 0.18 0.35 1.04 1.00 3526 2843
raceOtherHispanic 0.43 0.24 -0.07 0.89 1.00 4058 3114
raceNHBlack 0.54 0.15 0.24 0.82 1.00 3597 3177
raceOtherDMulti 0.82 0.19 0.45 1.19 1.00 3763 3257
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).p5 <- plot(bayes_fit) # Posterior distributions
p5 [[1]]
[[2]]
p6 <- pp_check(bayes_fit) # Posterior predictive checks
mcmc_trace(bayes_fit) # Convergence (optional)
bayes_R2(bayes_fit) # Model fit Estimate Est.Error Q2.5 Q97.5
R2 0.1313342 0.01265055 0.1064607 0.156078 # Leave-one-out cross-validation
ggsave("bayes_fit_plot.png", plot = p5, width = 8, height = 6, dpi = 300)
p6
ggsave("pp_check.png", plot = p6, width = 8, height = 6, dpi = 300)Below are reported odds ratio from the posterior predicted values and the Bayesian regression summary
# Posterior ORs (drop intercept, clean labels)
bayes_or <- posterior_summary(bayes_fit, pars = "^b_") %>%
as.data.frame() %>%
tibble::rownames_to_column("raw") %>%
dplyr::mutate(
term = gsub("^b_", "", raw),
term = gsub("race", "race:", term),
term = gsub("sex", "sex:", term),
term = gsub("OtherDMulti", "Other/Multi", term),
term = gsub("OtherHispanic", "Other Hispanic", term),
OR = exp(Estimate),
LCL = exp(Q2.5),
UCL = exp(Q97.5)
) %>%
dplyr::select(term, OR, LCL, UCL) %>%
dplyr::filter(term != "Intercept")
knitr::kable(
bayes_or %>%
dplyr::mutate(dplyr::across(c(OR,LCL,UCL), ~round(.x, 2))),
digits = 2,
caption = "Bayesian posterior odds ratios (95% CrI) — reference: NH White (race), Male (sex)"
)| term | OR | LCL | UCL |
|---|---|---|---|
| age_c | 2.99 | 2.64 | 3.37 |
| bmi_c | 1.87 | 1.71 | 2.05 |
| sex:Female | 0.52 | 0.42 | 0.63 |
| race:MexicanAmerican | 2.00 | 1.41 | 2.84 |
| race:Other Hispanic | 1.54 | 0.93 | 2.43 |
| race:NHBlack | 1.71 | 1.28 | 2.27 |
| race:Other/Multi | 2.27 | 1.56 | 3.28 |
# Combined table
if (!dir.exists("outputs")) dir.create("outputs", recursive = TRUE)
saveRDS(svy_fit, "outputs/svy_fit.rds")
saveRDS(pool_mi, "outputs/pool_mi.rds")
saveRDS(bayes_fit, "outputs/bayes_fit.rds")
saveRDS(svy_or, "outputs/survey_OR_table.rds")
saveRDS(mi_or, "outputs/mi_OR_table.rds")
saveRDS(bayes_or, "outputs/bayes_OR_table.rds")# Results
# ---- Build compact results table (BMI & Age only) ----
library(dplyr);
library(tidyr);
library(knitr);
library(stringr)
# pretty "OR (LCL–UCL)" string
fmt_or <- function(or, lcl, ucl, digits = 2) {
paste0(
formatC(or, format = "f", digits = digits), " (",
formatC(lcl, format = "f", digits = digits), "–",
formatC(ucl, format = "f", digits = digits), ")"
)
}
# guardrails: require these to exist from Modeling
stopifnot(exists("svy_or"), exists("mi_or"), exists("bayes_or"))
for (nm in c("svy_or","mi_or","bayes_or")) {
if (!all(c("term","OR","LCL","UCL") %in% names(get(nm)))) {
stop(nm, " must have columns: term, OR, LCL, UCL")
}
}
svy_tbl <- svy_or %>% mutate(Model = "Survey-weighted MLE")
mi_tbl <- mi_or %>% mutate(Model = "MICE pooled")
bayes_tbl <- bayes_or %>% mutate(Model = "Bayesian")
all_tbl <- bind_rows(svy_tbl, mi_tbl, bayes_tbl) %>%
mutate(term = case_when(
str_detect(term, "bmi_c|\\bBMI\\b") ~ "BMI (per 1 SD)",
str_detect(term, "age_c|\\bAge\\b") ~ "Age (per 1 SD)",
TRUE ~ term
)) %>%
filter(term %in% c("BMI (per 1 SD)", "Age (per 1 SD)")) %>%
mutate(OR_CI = fmt_or(OR, LCL, UCL, digits = 2)) %>%
select(Model, term, OR_CI) %>%
arrange(
factor(Model, levels = c("Survey-weighted MLE","MICE pooled","Bayesian")),
factor(term, levels = c("BMI (per 1 SD)","Age (per 1 SD)"))
)
res_wide <- all_tbl %>%
pivot_wider(names_from = term, values_from = OR_CI) %>%
rename(
`BMI (per 1 SD) OR (95% CI)` = `BMI (per 1 SD)`,
`Age (per 1 SD) OR (95% CI)` = `Age (per 1 SD)`
)
kable(
res_wide,
align = c("l","c","c"),
caption = "Odds ratios (per 1 SD) with 95% CIs across models"
)| Model | BMI (per 1 SD) OR (95% CI) | Age (per 1 SD) OR (95% CI) |
|---|---|---|
| Survey-weighted MLE | 1.89 (1.65–2.15) | 3.03 (2.70–3.40) |
| MICE pooled | 1.73 (1.58–1.89) | 2.90 (2.60–3.24) |
| Bayesian | 1.87 (1.71–2.05) | 2.99 (2.64–3.37) |
# Posterior predictive draws
#Posterior predictive checks (binary outcome)
pp_samples <- posterior_predict(bayes_fit, ndraws = 500) # 500 draws
# Check dimensions
dim(pp_samples) # rows = draws, cols = observations[1] 500 5592# Plot overlay of observed vs predicted counts (duplicate image)
ppc_dens_overlay(y = adult_imp1$diabetes_dx, yrep = pp_samples[1:50, ]) +
labs(title = "Posterior Predictive Check: Density Overlay") +
theme_minimal()
ppc_bars(y = adult_imp1$diabetes_dx, yrep = pp_samples[1:50, ])
# Alternative PP plots (histogram / barplot) for binary outcome (bar chart preferred being discrete outcome)
p7 <- ppc_bars(y = adult_imp1$diabetes_dx, yrep = pp_samples[1:50, ]) +
labs(title = "Posterior Predictive Check: Barplot of Counts") +
theme_minimal()
ggsave("ppc_bars.png", plot = p7, width = 8, height = 6, dpi = 300)
p7
#PP check for proportions (useful for binary) # mean comparison
## to check if the simulated means match the observed mean
## mean
p8 <- ppc_stat(y = adult_imp1$diabetes_dx, yrep = pp_samples[1:100, ], stat = "mean") +
labs(title = "Posterior Predictive Check: Mean of Replicates") +
theme_minimal()
ggsave("ppc_stat_mean.png", plot = p8, width = 8, height = 6, dpi = 300)
p8 
## sd
p9 <- ppc_stat(y = adult_imp1$diabetes_dx, yrep = pp_samples[1:100, ], stat = "sd") +
labs(title = "PPC: Standard Deviation of Replicates") +
theme_minimal()
ggsave("ppc_stat_sd.png", plot = p9, width = 8, height = 6, dpi = 300)
p9
# PP checks with bayesplot options
color_scheme_set("blue")
p10 <- ppc_scatter_avg(y = adult_imp1$diabetes_dx, yrep = pp_samples[1:100, ]) +
labs(title = "Observed vs Predicted (Avg) Posterior Predictive")
ggsave("ppc_scatter_avg.png", plot = p10, width = 8, height = 6, dpi = 300)
p10
library(brms)
library(dplyr)
# Posterior summary
post_sum <- posterior_summary(bayes_fit)
colnames(post_sum)[1] "Estimate" "Est.Error" "Q2.5" "Q97.5" library(posterior)
library(bayesplot)
# Extract posterior draws as a draws_df # simulate posterior outcomes
post <- as_draws_df(bayes_fit)
# Check parameter names
colnames(post) [1] "b_Intercept" "b_age_c" "b_bmi_c"
[4] "b_sexFemale" "b_raceMexicanAmerican" "b_raceOtherHispanic"
[7] "b_raceNHBlack" "b_raceOtherDMulti" "Intercept"
[10] "lprior" "lp__" ".chain"
[13] ".iteration" ".draw" # Density overlay for age and bmi
p11 <- mcmc_areas(post, pars = c( "b_age_c","b_bmi_c","b_sexFemale","b_raceMexicanAmerican", "b_raceOtherHispanic","b_raceNHBlack","b_raceOtherDMulti" ))
ggsave("mcmc_areas.png", plot = p11, width = 8, height = 6, dpi = 300)
p11
###
predicted <- fitted(bayes_fit, summary = TRUE)
observed <- adult_imp1[, c("bmi", "age")]
# Plot for **bmi** (obs vs pred)
library(ggplot2)
p12 <- ggplot(data = NULL, aes(x = observed$bmi, y = predicted[, "Estimate"])) +
geom_point() +
geom_errorbar(aes(ymin = predicted[, "Q2.5"], ymax = predicted[, "Q97.5"])) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
xlab("Observed bmi") + ylab("Predicted bmi")
ggsave("plot_bmi_obs_pred.png", plot = p12, width = 8, height = 6, dpi = 300)
p12
# Combine observed and predicted into one data frame
plot_data <- adult_imp1 %>%
mutate(
predicted_bmi = predicted[, "Estimate"],
lower_ci = predicted[, "Q2.5"],
upper_ci = predicted[, "Q97.5"],
obs_index = 1:nrow(adult_imp1) # index for x-axis
)
# Line plot
p13 <- ggplot(plot_data, aes(x = obs_index)) +
geom_line(aes(y = bmi, color = "Observed")) + # observed BMI
geom_line(aes(y = predicted_bmi, color = "Predicted")) + # predicted BMI
geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.2) + # uncertainty
labs(x = "Observation", y = "BMI", color = "Legend") +
theme_minimal()
ggsave("Line plot_bmi_obs_pred.png", plot = p13, width = 8, height = 6, dpi = 300)
p13
###
summary(adult_imp1$bmi) Min. 1st Qu. Median Mean 3rd Qu. Max.
14.1 24.1 27.7 29.0 32.4 82.9 summary(plot_data$bmi_c) Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.09476 -0.69669 -0.19338 -0.01133 0.46371 7.52398 prior_summary(bayes_fit) prior class coef group resp dpar nlpar lb ub
normal(0, 2.5) b
normal(0, 2.5) b age_c
normal(0, 2.5) b bmi_c
normal(0, 2.5) b raceMexicanAmerican
normal(0, 2.5) b raceNHBlack
normal(0, 2.5) b raceOtherDMulti
normal(0, 2.5) b raceOtherHispanic
normal(0, 2.5) b sexFemale
student_t(3, 0, 10) Intercept
tag source
user
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
userprior_draws <- tibble(
term = rep(c("BMI (per 1 SD)", "Age (per 1 SD)"), each = 4000),
estimate = c(rnorm(4000, 0, 1), rnorm(4000, 0, 1)),
type = "Prior"
)
summary(prior_draws) term estimate type
Length:8000 Min. :-3.585220 Length:8000
Class :character 1st Qu.:-0.696978 Class :character
Mode :character Median : 0.000572 Mode :character
Mean :-0.007125
3rd Qu.: 0.668627
Max. : 3.922817 post# A draws_df: 1000 iterations, 4 chains, and 11 variables
b_Intercept b_age_c b_bmi_c b_sexFemale b_raceMexicanAmerican
1 -2.6 1.1 0.70 -0.71 0.67
2 -2.7 1.0 0.62 -0.57 0.65
3 -2.6 1.1 0.65 -0.76 0.63
4 -2.7 1.0 0.65 -0.67 0.82
5 -2.6 1.1 0.61 -0.73 0.75
6 -2.5 1.0 0.60 -0.77 0.61
7 -2.8 1.1 0.66 -0.66 0.52
8 -2.8 1.2 0.67 -0.57 0.94
9 -2.8 1.1 0.65 -0.52 0.84
10 -2.6 1.1 0.67 -0.85 0.70
b_raceOtherHispanic b_raceNHBlack b_raceOtherDMulti
1 0.605 0.52 0.95
2 0.338 0.45 0.69
3 0.566 0.63 0.54
4 0.453 0.61 0.78
5 0.090 0.50 0.62
6 0.015 0.48 0.60
7 0.736 0.50 0.84
8 0.913 0.57 1.07
9 0.570 0.66 0.81
10 0.467 0.54 0.97
# ... with 3990 more draws, and 3 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}head(post)# A draws_df: 6 iterations, 1 chains, and 11 variables
b_Intercept b_age_c b_bmi_c b_sexFemale b_raceMexicanAmerican
1 -2.6 1.1 0.70 -0.71 0.67
2 -2.7 1.0 0.62 -0.57 0.65
3 -2.6 1.1 0.65 -0.76 0.63
4 -2.7 1.0 0.65 -0.67 0.82
5 -2.6 1.1 0.61 -0.73 0.75
6 -2.5 1.0 0.60 -0.77 0.61
b_raceOtherHispanic b_raceNHBlack b_raceOtherDMulti
1 0.605 0.52 0.95
2 0.338 0.45 0.69
3 0.566 0.63 0.54
4 0.453 0.61 0.78
5 0.090 0.50 0.62
6 0.015 0.48 0.60
# ... with 3 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}names(prior_draws)[1] "term" "estimate" "type" library(brms)
library(tidyr)
# Extract posterior draws
post <- as_draws_df(bayes_fit) %>% # bayes_fit = your brms model
select(b_bmi_c, b_age_c) %>% # select your coefficient columns
pivot_longer(
everything(),
names_to = "term",
values_to = "estimate"
) %>%
mutate(
term = case_when(
term == "b_bmi_c" ~ "BMI (per 1 SD)",
term == "b_age_c" ~ "Age (per 1 SD)"
),
type = "Posterior"
)
## visualization of prior and predicted draws
combined_draws <- bind_rows(prior_draws, post)
library(ggplot2)
p14 <- ggplot(combined_draws, aes(x = estimate, fill = type)) +
geom_density(alpha = 0.4) +
facet_wrap(~ term, scales = "free", ncol = 2) +
theme_minimal(base_size = 13) +
labs(
title = "Prior vs Posterior Distributions",
x = "Coefficient estimate",
y = "Density",
fill = ""
)
p14
ggsave("Prior vs Posterior Distributions_bmi_age.png", plot = p14, width = 8, height = 6, dpi = 300)
#### Compute proportion of diabetes=1 for each draw
pp_proportion <- rowMeans(pp_samples) # proportion of 1's in each posterior draw
# Summary of posterior proportions
summary(pp_proportion) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.09371 0.10479 0.10908 0.10923 0.11373 0.13215 # Optional: visualize the posterior probability distribution
pp_proportion_df <- tibble(proportion = pp_proportion)
p15 <- ggplot(pp_proportion_df, aes(x = proportion)) +
geom_histogram(binwidth = 0.01, fill = "skyblue", color = "black") +
labs(
title = "Posterior Distribution of Proportion of Diabetes = 1",
x = "Proportion of Diabetes = 1",
y = "Frequency"
) +
theme_minimal()
p15
ggsave("Posterior Distribution of Proportion of Diabetes = 1.png", plot = p15, width = 8, height = 6, dpi = 300)####
library(tidyverse)
# Posterior predicted proportion vector
# pp_proportion <- rowMeans(pp_samples) # if not already done
known_prev <- 0.089 # NHANES prevalence
# Posterior summary
posterior_mean <- mean(pp_proportion)
posterior_ci <- quantile(pp_proportion, c(0.025, 0.975)) # 95% credible interval
# Create a data frame for plotting
pp_df <- tibble(proportion = pp_proportion)
# Plot
p16 <- ggplot(pp_df, aes(x = proportion)) +
geom_histogram(binwidth = 0.005, fill = "skyblue", color = "black") +
geom_vline(xintercept = known_prev, color = "red", linetype = "dashed", size = 1) +
geom_vline(xintercept = posterior_mean, color = "blue", linetype = "solid", size = 1) +
geom_rect(aes(xmin = posterior_ci[1], xmax = posterior_ci[2], ymin = 0, ymax = Inf),
fill = "blue", alpha = 0.1, inherit.aes = FALSE) +
labs(
title = "Posterior Predicted Diabetes Proportion vs NHANES Prevalence",
subtitle = paste0("Red dashed = NHANES prevalence (", known_prev,
"), Blue solid = Posterior mean (", round(posterior_mean,3), ")"),
x = "Proportion of Diabetes = 1",
y = "Frequency"
) +
theme_minimal()
ggsave("Posterior Distribution of Proportion of Diabetes = 1.png", plot = p16, width = 8, height = 6, dpi = 300)
p16
library(dplyr)
# Posterior predicted proportion
posterior_mean <- mean(pp_proportion)
posterior_ci <- quantile(pp_proportion, c(0.025, 0.975)) # 95% credible interval
# NHANES prevalence with SE from survey::svymean
# Suppose you already have:
# svymean(~diabetes_dx, nhanes_design_adult, na.rm = TRUE)
known_prev <- 0.089 # Mean prevalence
known_se <- 0.0048 # Standard error from survey
# Calculate 95% confidence interval
known_ci <- c(
known_prev - 1.96 * known_se,
known_prev + 1.96 * known_se
)
# Print results
data.frame(
Type = c("Posterior Prediction", "NHANES Prevalence"),
Mean = c(posterior_mean, known_prev),
Lower_95 = c(posterior_ci[1], known_ci[1]),
Upper_95 = c(posterior_ci[2], known_ci[2])
) Type Mean Lower_95 Upper_95
2.5% Posterior Prediction 0.109235 0.0975456 0.1207082
NHANES Prevalence 0.089000 0.0795920 0.0984080####
library(ggplot2)
library(dplyr)
# Create a data frame for plotting
ci_df <- data.frame(
Type = c("Posterior Prediction", "NHANES Prevalence"),
Mean = c(0.1096674, 0.089),
Lower_95 = c(0.09772443, 0.079592),
Upper_95 = c(0.1210658, 0.098408)
)
# Plot
p17 <- ggplot(ci_df, aes(x = Type, y = Mean, color = Type)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = Lower_95, ymax = Upper_95), width = 0.2) +
ylim(0, max(ci_df$Upper_95) + 0.02) +
labs(
title = "Comparison of Posterior Predicted Diabetes Proportion vs NHANES Prevalence",
y = "Proportion of Diabetes",
x = ""
) +
theme_minimal(base_size = 14) +
theme(legend.position = "none")
ggsave("Population Prevalence vs Proportion of Diabetes = 1.png", plot = p17, width = 8, height = 6, dpi = 300)
p17
Multiple Linear Regression
Model: svyglm(formula = form_cc, design = des_cc, family = quasibinomial())
All predictors are significant: age (p < 0.001) and BMI (p < 0.001) show strong positive associations with the outcome, while being female (p = 0.0004) is negatively associated. Other significant associations include raceMexican American (p = 0.0008), raceOther Hispanic (p = 0.0087), raceNH Black (p = 0.0117), and raceOther/Multi (p = 0.0014).
MICE
All predictors are statistically significant.
Positive associations: age (p < 0.001), BMI (p < 0.001), and all race categories compared to reference.
Negative association: being female (p < 0.001)
Bayesian Regression
Sampling: NUTS (4 chains, 2000 iterations each; 1000 warmup, 4000 post-warmup draws)
Convergence & Diagnostics
Rhat = 1.00 for all parameters → excellent convergence
Bulk_ESS / Tail_ESS: Large values (>2000) → high effective sample sizes, reliable posterior estimates.
Interpretation
Strong predictors: Age and BMI are strongly positively associated with diabetes risk.
Sex effect: Females have a lower probability of diabetes compared to males
R² = 0.13 shows 13% of the variance in the outcome (diabetes_dx) is explained by your predictors (age, BMI, sex, race), 95% Credible Interval: 0.106–0.156, indicates that, given your model and data, the true proportion of variance explained is plausibly between 10.6% and 15.6% and shows uncertainty in the explained variance, which is natural for probabilistic models.
Est.Error = 0.0127, reflects the standard error of the R² estimate across posterior samples. The small SE indicates that the R² estimate is fairly precise.
age_c-Each 1-unit increase in centered age increases the log-odds of diabetes by 1.09. Strong positive effect.
bmi_c-Higher BMI is associated with higher diabetes risk. sexFemale-Females have lower odds of diabetes compared to males.
raceMexicanAmerican-Higher odds of diabetes vs. reference race (likely NH White)
raceOtherHispanic-Slightly higher odds vs reference, but interval crosses zero → uncertain effect.
raceNHBlack-Significantly higher odds of diabetes compared to reference. raceOtherDMulti-Higher odds of diabetes vs reference group.
Posterior distribution of all parameters in the model. (1)Density plot of posterior samples each parameter (e.g., intercept, slope) into a smoothed density curve, showing most of the posterior probability mass lies for best estimates and uncertainty.
Posterior Predictive Distribution - generated from posterior predictive draws: 𝑦rep∼𝑝(𝑦new∣𝜃)yrep∼p(ynew∣θ) simulate the data given posterior parameter estimates.Posterior predictive checks (PPC) compare these simulations to real data to assess model fit.
Incorporating Uncertainty two sources of uncertainty: Parameter uncertainty: captured in the posterior distributions Predictive uncertainty: captured in posterior predictive draws
Combining the two provide credible intervals for predictions, not just point predictions and specifies - Given the BMI, the probability of diabetes is likely between 40–55%.”
Comparing Models
All three models (survey-weighted MLE, multiple imputation, Bayesian) agree closely on the direction and magnitude of the effects of BMI and age.
Age is a stronger predictor than BMI, nearly tripling the odds per 1 SD.
BMI significantly increases diabetes risk (~1.7–1.9× per 1 SD).
Differences between models are minor, indicating robust and reliable findings despite missing data or modeling approach.
The difference could be due to:
The results find our model is reasonable but slightly conservative (predicting higher risk) relative to the observed population prevalence.
Across multiple modeling approaches—survey-weighted maximum likelihood, multiple imputation, and Bayesian regression—both age and BMI were consistently strong predictors of diabetes. Eachstandard deviation increase in age nearly tripled the odds of diabetes, while a similar increase in BMI elevated the odds by approximately 1.7–1.9 times. The consistency of these results across models highlights the robustness of the associations and underscores the importance of age and BMI as key risk factors for diabetes in this population.
Effect stability: point estimates in rhe Bayesian model’s closely aligned with those from the frequentist, indicating that prior regularization stabilized the estimates in the presence of modest missingness.
Uncertainty quantification: Bayesian credible intervals of odds ration were slightly narrower yet overlapped the frequentist confidence intervals, suggest comparable inferential precision while offering improved interpretability.
Design considerations: # Survey-weighted MLE (Maximum Likelihood Estimator) - incorporates each observation weighted according to its survey weight. - provide estimates that reflect the population-level parameters, not just the sample- produces population-representative estimates. # Bayesian model with normalized weights- - instead of fully modeling the survey design, it used normalized sampling weights as importance weights - the scaled weights that sum to the sample size approximates the effect of survey weights, but does not fully account for: Stratification, clustering, design-based variance adjustments. - Bayesian inference treats the weighted likelihood as from a regular model, ignoring some survey design features.
The use of multiple imputation allowed for robust analysis despite missing data, increasing power and reducing bias. Comparison of frequentist and Bayesian models demonstrated consistency in significant predictors, while Bayesian approaches provided the advantage of posterior distributions and probabilistic interpretation. The =
Across all models, both age and BMI emerged as strong and consistent predictors of diabetes. The consistency across modeling approaches strengthens the validity of these findings Multiple imputation accounted for potential biases due to missing data, and Bayesian modeling provided robust credible intervals that closely matched frequentist estimates. align with previous epidemiological research indicating that increasing age and higher BMI are among the most important determinants of type 2 diabetes risk.Cumulative exposure to metabolic and lifestyle risk factors over time, and the role of excess adiposity and insulin related effects account for diabetes.
Survey weighted dataset strenghthens ensuring population representativeness, multiple imputation to handle missing data, and rigorous Bayesian estimation provided high effective sample sizes and R̂ ≈ 1.00 across parameters confirmed excellent model convergence. Bayesian logistic regression provided inference statistically consistent and interpretable achieving the aim of this study. In future research hierarchical model using NHANES cycles or adding variables (lab tests) could assess nonlinear effects of metabolic risk factors.
Our study was a cross-sectional study design - precludes potential residual confounding from unmeasured factors such as diet, physical activity, and genetic predisposition